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Abstract 

We present a novel and powerful Compressible High-ORder Unstructured Spectral-difference (CHO¬ 
RUS) code for simulating thermal convection and related fluid dynamics in the interiors of stars 
and planets. The computational geometries are treated as rotating spherical shells filled with strat¬ 
ified gas. The hydrodynamic equations are discretized by a robust and efficient high-order Spectral 
Difference Method (SDM) on unstructured meshes. The computational stencil of the spectral 
difference method is compact and advantageous for parallel processing. CHORUS demonstrates 
excellent parallel performance for all test cases reported in this paper, scaling up to 12,000 cores on 
the Yellowstone High-Performance Computing cluster at NCAR. The code is verified by defining 
two benchmark cases for global convection in Jupiter and the Sun. CHORUS results are compared 
with results from the ASH code and good agreement is found. The CHORUS code creates new op¬ 
portunities for simulating such varied phenomena as multi-scale solar convection, core convection, 
and convection in rapidly-rotating, oblate stars. 

Keywords: Spectral difference method, High-order, Unstructured grid, Astrophysical fluid 
dynamics 


1. Introduction 

Turbulent convection is ubiquitous in stars and planets. In intermediate-mass stars like the Sun, 
convection acts with radiation to transport the energy generated by fusion in the core to the surface 
where it is radiated into space. In short, convection enables the Sun to shine. It also redistributes 
angular momentum, establishing differential rotation (equator spinning about 30% faster than the 
polar regions) and meridional circulations (with poleward flow near the surface). Furthermore, 
turbulent solar convection and the mean flows it establishes act to amplify and organize magnetic 
fields, giving rise to patterns of magnetic activity such as the 11-year sunspot cycle. Other stars 
similarly exhibit magnetic activity that is highly correlated with the presence of surface convection 
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and differential rotation minis]. Stars are hydromagnetic dynamos, generating vibrant, sometimes 
cyclic, magnetic activity from the kinetic energy of plasma motions. 

Perhaps the biggest challenge in modeling solar and stellar convection is the vast range of spatial 
and temporal scales involved. Solar observations reveal a network of convection cells on the surface 
of the Sun known as granulation [3]. Each cell has a characteristic size of about 1000 km and a 
lifetime of 10-15 min. However, in order to account for the differential rotation and cyclic magnetic 
activity of the Sun, larger-scale convective motions must also be present, occupying the bulk of the 
convection zone that extends from the surface down to 0.7 Rq where Rq is the solar radius [5]. 
These so-called “giant cells” have characteristic length and time scales of order 100,000 and several 
weeks respectively. Below the convective zone lies the convectively stable radiative zone where 
energy is transported by radiative diffusion. The interface between the convective and radiative 
zones is a thin internal boundary layer that poses its own modeling challenges. Here convection 
overshoots into the stable interior, exciting internal gravity waves and establishing a layer of strong 
radial shear in the differential rotation that is known as the solar tachocline [5]. 

Another formidable modeling challenge is the geometry. High-mass stars (M > 4M 0 , where 
Mq is the solar mass) are inverted Suns, with convectively stable (radiative) envelopes surrounding 
convectively unstable cores. These are also expected to possess vigorous dynamo action but much 
of it is likely hidden from us, occurring deep below the surface where it cannot be observed with 
present techniques. Stellar lifetimes are anticorrelated with mass, so all high-mass stars are signifi¬ 
cantly younger than the Sun. Furthermore, stars spin down as they age due to torques exerted by 
magnetized stellar winds mmm- Thus, most high-mass stars spin much faster that the Sun, as 
much as one to two orders of magnitude. The fastest rotators are significantly oblate. For example, 
the star Regulus in the constellation of Leo has an equatorial diameter that is more than 30% larger 
than its polar diameter {9j. 

Other types of stars and planets pose their own set of challenges. Low-mass main sequence stars 
(M < M 0 ) are convective throughout, from their core to their surface, so they require modeling 
strategies that can gracefully handle the coordinate singularity at the origin (r = 0) as well as the 
the small-scale convection established by the steep density stratification and strong radiative cooling 
in the surface layers. Red giants have deep convective envelopes and dense, rapidly-rotating cores. 
Jovian planets have both deep convective dynamos and shallow, electrically neutral atmospheric 
dynamics that drive strong zonal winds. 

All of these systems have the common property that they are highly turbulent. In other words, 
the Reynolds number R e = UL/v is very large, where U and L are velocity and length scales and 
v is the kinematic viscosity of the plasma. For example, in the solar convection zone R e > 10 12 
[5j. Furthermore, though all of these systems are strongly stratified in density (compressible), most 
possess convective motions that are much slower than the sound speed c s . Thus, the Mach number 
M a « 1. Exceptions include surface convection in solar-like stars and low-mass stars and deep 
convection in the relatively cool red giants, where M a can approach unity. 

Meeting this considerable list of challenges requires a flexible, accurate, robust, and efficient 
computational algorithm. It is within this context that we here introduce the Compressible High- 
Order Unstructured Spectral difference (CHORUS) code. CHORUS is the first numerical model of 
global solar, stellar, and planetary convection that uses an unstructured grid. This valuable feature 
allows CHORUS to avoid the spherical coordinate singularities at the origin (r = 0) and at the 
poles (colatitude 6 = 0, 7r) that plague codes based on structured grids in spherical coordinates. 
Such coordinate singularities compromise computational efficiency as well as accuracy, due to the 
grid convergence that can place a disproportionate number of points near the singularities and 
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that can severely limit the time step through the Courant-Freidrichs-Lewy (CFL) condition. Thus, 
CHORUS can handle a wide range of global geometries, from the convective envelopes of solar-like 
stars and red giants to the convective cores of high-mass stars. 

The flexibility of the unstructured grid promotes maximal computational efficiency for capturing 
multi-scale nature of solar and stellar convection. As noted above, boundary layers play an essential 
role in the internal dynamics of stars and planets. In solar-like stars with convective envelopes, 
much of the convective driving occurs in the surface layers, producing granulation that transitions 
to giant cells through a hierarchical merging of downflow plumes [5]. Meanwhile, the tachocline 
and overshoot region at the base of the convection zone play a crucial role in the solar dynamo. 
The unstructured grid of CHORUS will enable us to locally enhance the resolution in these regions 
in order to capture the essential dynamics. Similarly, optimal placing of grid points will allow 
CHORUS to efficiently model other phenomena such as core-envelope coupling in red giants (see 
above). Furthermore, the unstructured grid is deformable. So, it can handle the oblateness of 
rapidly-rotating stars as well as other steady and time-dependent distortions arising from radial 
pulsation modes or tidal forcing by stellar or planetary companions. 

The spectral difference method (SDM) employed for CHORUS achieves high accuracy for a 
wide range of spatial scales, which is necessary to capture the highly turbulent nature of stellar and 
planetary convection. Furthermore, its low intrinsic numerical dissipation enables it to be run in an 
inviscid mode (no explicit viscosity), thus maximizing the effective R e for a given spatial resolution. 

One feature of CHORUS that is not optimal for modeling stars is the fully compressible nature 
of the governing equations. The dynamics of stellar and planetary interiors typically operates at 
low Mach number such that acoustic time scales are orders of magnitude smaller than the time 
scales for convection or other large-scale instabilities. Therefore, a fully compressible solver such as 
CHORUS is limited by the CFL constraint imposed by acoustic waves. Many codes circumvent this 
problem by adopting anelastic or pseudo-compressible approximations that filter out sound waves. 
We choose instead to define idealized problems by scaling up the luminosity to achieve higher Mach 
numbers while leaving other important dynamical measures such as the Rossby number unchanged. 
This allows us to take advantage of the hyperbolic nature of the compressible equations, which is well 
suited for the SDM method and which promotes excellent scalability on massively parallel computing 
platforms. Furthermore, the compressible nature of the equations will enable CHORUS to address 
problems that are inaccessible or at best challenging with anelastic and pseudo-compressible codes. 
These include high-Mach number convection in Red Giants, the coupling of photospheric and deep 
convection, and the excitation of the radial and non-radial acoustic oscillations (p-modes) that form 
the basis of helio- and asteroseismology. 

Currently CHORUS employs an explicit time-stepping method, which is not optimal for low 
Mach number flow, especially on non-uniform grids. However, the SDM is well suited for split¬ 
timestepping and implicit time-stepping methods which we intend to implement in the future. 

The purpose of this paper is to introduce the CHORUS code, to describe its numerical algorithm, 
and to verify it by comparing it to the well-established Anelastic Spherical Harmonic (ASH) code. 
We begin with a discussion of the mathematical formulation and numerical algorithms of CHORUS 
in sections 2 and 3 respectively. In section 4 we generalize the anelastic benchmark simulations of 
m to provide initial conditions for compressible models like CHORUS. We then address boundary 
conditions in section 5 and the conservation of angular momentum in section 6, which can be a 
challenge for global convection codes. In section 7 we verify the CHORUS code by comparing its 
output to analogous ASH simulations for two illustrative benchmark cases, representative of Jovian 
planets and solar-like stars. We summarize our results and future plans in section 8. 
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2. Mathematical Formulation 


We consider a spherical shell of ideal gas, bounded by an inner spherical surface at r = Ri 
and an outer surface at r = R 0 where r is the radius. We assume that the bulk of the mass is 
concentrated between both surfaces, and the gravity satisfies g = —gr = — where G is the 
gravitational constant, M is the interior mass and r is the radial unit vector. Consider a reference 
frame that is uniformly rotating about the z axis with angular speed f2 0 = fl 0 e z where e z is the unit 
vector in z direction. In this rotating frame, the effect of Coriolis force is added to the momentum 
conservation equations. The Centrifugal forces are negligible as they have much less contribution in 
comparison with the gravity. However, the CHORUS code can handle oblate spheroids in rapidly- 
rotating objects and we will consider Centrifugal force effect in the future. The resulting system of 
hydrodynamic equations is 

! = -V.(„u), (1) 

= -V • puu Vp T V • r T pg 2pfl 0 x u, (2) 

at 

dE 

— = -V • ((E + p) u) + V • (u • t - f) + pu ■ g, (3) 

where f, p, T, p and u are time, pressure, temperature, density and velocity vector respectively. 
E is the total energy per unit volume and is defined as E = + \pvL ■ u where 7 is the ra¬ 

tio of the specific heats, r is the viscous stress tensor for a Newtonian fluid. The term u • r 
in Eq. © represents the viscous heating. The diffusive flux f is generally treated in the form of 
f = — KpT\7S — K r pC p \7T where k is the entropy diffusion coefficient, S is the specific entropy, n r 
is the thermal diffusivity (thermal conductivity and radiative conductivity), and C p is the specific 
heat at constant pressure. The entropy diffusion is a popular way of parameterizing the energy flux 
due to unresolved, subgrid-scale convective motions which tend to mix entropy Dunging. At the 
bottom of the convection zone, f is generally prescribed with a constant value which acts as the 
energy source of convection. The last term in Eq. is the work done by buoyancy. 


The governing equations for fully compressible model can be written in a conservative form as 


dQ <9F dG 
dt dx dy 


<9H 

dz 


+ M = 0, 


( 4 ) 


where Q is the vector of conserved variables, M is the combination of the Coriolis force term and 
the gravitational force term, and F, G, H are the total fluxes including both inviscid and viscous 
flux vectors in local Cartesian coordinates (transforming to an arbitrary geometry will be discussed 
in next section). We write these as F = F inv — F„, G = G i nv — G„, and H = H inv — H„, where 
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In Eq.([|)-([ 8 |, u, v, and w are the velocity components in the x, y , and z directions respectively. 
The viscous stress tensor components in Eq.([7| and ([ 8 ]) can be written in the following form 
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( 9 ) 


/ ^x i u y i \ 

t zz = 2p(w z ---), 


where p = pv is the dynamic viscosity and v is the kinematic viscosity. 


3. Numerical Algorithm 

The equations of motion (4)-(9) are solved using a Spectral Difference method (SDM). The SDM 
is similar to the multi-domain staggered method originally proposed by Kopriva and his colleagues 
[ 151 . Liu et al [2] first formulated the SDM for wave equations by extending the multi-domain 
staggered method to triangular elements. The SDM was then employed by Wang et al for inviscid 
compressible Euler equations on simplex elements [15] and viscous compressible flow on unstructured 
grids [16]. The SDM is simpler than the traditional Discontinuous Galerkin (DG) method [TTj since 
DG deals with the weak form of the equations and involves volume and surface integrals. The 
weak form is developed by integrating the product of a test function with the compressible Navier- 
Stokes equations. For the discontinuous Galerkin method, the integration is performed over the 
spatial coordinates of each finite element of the computational domain. The SDM is similar to 
the quadrature-free nodal discontinuous Galerkin method [181. The SDM of CHORUS is designed 
for unstructured meshes of all hexaliedral elements m ■ This spectral difference approach employs 
high-order polynomials within each hexahedral element locally. In particular, we employ the roots 
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(a) Full spherical shell mesh 


(b) Spherical shell segment spanning 1/8 of the volume 


Figure 1: Unstructured mesh consisting of all hexahedral elements. The full spherical shell mesh (a) is generated by 
fusing eight segments as illustrated in (b) using the GAMBIT software m- At present, nodal points are uniformly 
distributed in the radial direction but GAMBIT also allows for variable radial resolution. 


of Legendre polynomials plus two end points for locating flux points |20| . The stability for linear 
advection of this type of SDM has been proven by Jameson [2T]. Overall, high-order SDM is simple 
to formulate and CHORUS is very suitable for massively parallel processing. 

The computational domain is divided into a collection of non-overlapping hexahedral elements 
as illustrated in Fig [I] These elements share similarities with the control volumes in the Finite 
Volume Method. To achieve an efficient implementation, all hexahedral elements in the physical 
domain (x,y,z) are transformed to a standard cube (0 < £ < 1, 0 < 77 < 1, 0 < ( < 1). This 
mapping is achieved through a Jacobian matrix 


d(x,y,z) 


X£ Xyj Xq 

Vi Vv VC 

Z£ z v 


( 10 ) 


The governing equations in conservative form in the physical domain as described by Eq.Q are 
then transformed into the computational domain. The transformed equations are written as 


dQ dF dG <9H 

— 1 -I-1-1-h 

dt d£ dr] d( 


M = 0, 


( 11 ) 


where Q = | J7|Q and M = | J'lM using the determinant \ J\ of J. The transformed flux components 
can be written as a combination of the physical flux components as 


F 
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( 12 ) 
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In the two-dimensional standard element as illustrated in Fig(2j two sets of points are defined, 
namely the solution and flux points for the SDM. A total of 9 solution points and 24 flux points are 
employed for the third-order SDM in 2D. A more detailed description of the SDM for quadrilateral 
elements can be found in [[23] . The unknown conserved variables are stored at the solution points, 
while flux components F, G, and H are stored at the flux points in corresponding directions. 
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Flux points 9 Solution points 

Figure 2: Layout of a standard element for a 3rd order SDM in 2D 


In order to construct a degree (N-l) polynomial in each coordinate direction, N solution points 
are required (thus N is defined as the order the scheme). The solution points in ID are chosen to 
be Chebyshev-Gauss-quadrature points defined by 


X B 
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s = 1,2, - -- ,N. 


(13) 


The flux points are selected to be the Legendre-Gauss-quadrature points plus the two end 
points, 0 and 1. Choosing P-i(£) = 0 and Po(£) = 1, we can determine the higher-degree Legendre 
polynomials as 

O77 — 1 77 — 1 

Pn{0 = — -(2£-l)P„_i(0--- Pn- 2 ( 0 . n =!,■■■ ,N —1. (14) 

n n 

The locations of these Legendre-Gauss quadrature points for the N-th order SDM are the roots 
of equation P n _i(£) plus two end points. 

Using the solutions at N solution points, a degree (N-l) polynomial can be built using the 
following Lagrange basis defined as 


hi(X) 
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Similarly, using the (N+l) fluxes at the flux points, a degree N polynomial can be built for the 
flux using a similar Lagrange basis defined as 


h + jW 
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( 16 ) 


The reconstructed solution for the conserved variables in the standard element is just the tensor 
products of the three one-dimensional polynomials, i.e., 


N N N jz 
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Similarly, the reconstructed flux polynomials take the following forms: 
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The flux polynomials are element-wise continuous, but discontinuous across element interfaces. 
For computing the inviscid fluxes, an approximate Riemann solver l24tl25| is employed to compute a 
common flux at interfaces and to ensure conservation and stability. Here the Rusanov flux treatment 
for the £ direction is formulated as F rus = 1(11'™’' + F 1 ™ — sgn(n ■ V£)(H4| + c s )(Qr — Ql)\ </V£|), 
where n is the interface normal direction, V n is the fluid velocity normal to the interface and c s 
is the speed of sound. If the normal direction of the cell interface is mapped to either the 77 or £ 
direction, the Riemann fluxes can be formulated similarly. 

For calculating the viscous fluxes, a simple averaging procedure is used for evaluating fluxes at 
interfaces [23]. This procedure is similar to the BRI scheme [26]. For future implementation of 
implicit time stepping methods, we can extend the CHORUS code to use BR2 scheme m- 

The number of Degrees of Freedom (DOFs) for CHORUS simulations in this paper is computed 
as 

DOFs = N e i ement x N 3 , (19) 

where N e i ement is the total number of elements in the spherical shell and N is the order of the 
scheme, which is equal to the number of solution points in each direction within one standard 
element. 

The CHORUS code is written in FORTRAN 90 and efficient parallel performance is achieved by 
using the Message Passing Interface (MPI) for interprocessor communication. The ParMetis package 
[25] is used to partition the unstructured mesh by means of a graph partitioning method. The 
parallel scalability of the CHORUS code is shown in Fig[3]using the Yellowstone High-Performance 
computing cluster at the National Center for Atmospheric Research (NCAR), which is built on 
IBM’s iDataPlex architecture with Intel Sandy Bridge processors. These numerical experiments 




demonstrate strong scaling, with Tf denoting the execution time of the sequential CHORUS code 
and T n denoting the execution time of the parallel CHORUS code with n processors. Two sets 
of simulations using the 4th order SDM are shown. The total numbers of elements in physical 
domain for testl and test2 are 294,912 and 1,105,920 respectively, which correspond to 18,874,368 
and 70,778,880 DOFs in the 4th-order SDM simulations. 



Figure 3: Strong scaling of the CHORUS code on Yellowstone for 18.9M (testl) and 70.8M (test2) DOFs. The latter 
test in particular achieves over 92% effeciency for 12k processors. 


4. Initial Conditions and Stratification 

Though the equations that CHORUS solves are general, we are intersted in simulating global- 
scale convection in stellar and planetary interiors for reasons discussed in Section 1. Thermal 
convection is a classical fluid instability in the sense that it can develop from a static equilibrium 
state that satisfies certain instability criteria [29] ■ The first is the Schwarzschild criterion which 
requires a negative (superadiabatic) radial entropy gradient dS/dr < 0. The second is that the 
buoyancy force must be sufficiently strong to overcome viscous and thermal diffusion. This is 
typically quantified in terms of the Rayleigh number, which must exceed a critical value in order 
for convection to ensue. 

Though linear theory is generally concerned with static, equilibrium initial conditions, numerical 
simulations can tolerate initial conditions that are not in equilibrium. However, it is still desirable 
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to initiate nonlinear simulations with states that are close to equilibrium in order to mitigate violent 
initial transients and minimize nonlinear (dynamic) equilibration times. 

In the sections that follow we describe how we set up the initial conditions for spherical shells 
of convection. CHORUS can also handle convective cores and fully convective geometries but we 
defer these applications to future papers. Note that these initial conditions are static relative to 
the uniform rotation of the coordinate system (u = 0). Thus, the differential rotation, meridional 
circulation, and convective motions are in no way imposed; they are zero initially. After specifying 
the static initial conditions, CHORUS automatically introduces random thermal perturbations 
through the non-axisymmetric distribution of unstructured grid points to excite the convection 
which in turn establishes the mean flows. Note also that the stratification in a stellar or planetary 
convection zone is nearly hydrostatic. Convection will modify this but only slightly. So, the initial 
conditions not only excite the convection but they also establish the basic background stratification 
including crucial simulation properties such as the density contrast across the convection zone. 


4-1. Static Equilibrium Equations 

To establish the static initial conditions we first consider a steady state in the absence of motions 
( d/dt = 0, u = 0). Then the governing equations 0-0 reduce to the equation of hydrostatic 
balance, 


dpo 

dr 

and the equation of thermal energy balance, 


~Po9, 


( 20 ) 


d 

dr 


(' r 2 K.p 0 T 0 


dSo 

dr 


r K r 


dJ o 
dr 


) = o. 


( 21 ) 


where the subscript 0 denotes the initial state. These are supplemented with the ideal gas law 


Po = WpoTo, 

where IK is the specific gas constant, and the equation for specific entropy, 

1/7 

S 0 =C p \n(^). 

P 0 


( 22 ) 


(23) 


4-2. Polytropic, Adiabatic Reference State 

Though we solve the equations in dimensional form, it is useful to define several nondimensional 
numbers that characterize the parameter regime of the solution: 


_ GMdAS v 

^ 1 r — 5 ■E'k — 

VUC n K 


v at 1 f Pi\ a Ri 

——,AT p = M —),P = ~b~- 

p 0 Rq 


(24) 


Here R a is the Rayleigh number, P r is the fluid Prandtl number, Ek is the Ekman number, exp(iVp) 
is the density ratio across the layer, with pi and p c as the the densities at the inner and outer 
boundaries, d = R 0 — Ri, /3 is the aspect ratio, and AS is the entropy difference across the layer, 
averaged over latitude and longitude. 

The hydrostatic balance equation (20), along with the constitutive equations ( f22j ) and (23), can 
be satisfied by introducing a polytropic stratification as described by Jones et al. [10]: 

Pa = PcX n , T a = T c \, and p a = p c x n+1 , (25) 
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where 


(26) 
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(27) 


The subscripts i, o , and c refer to the bottom, top, and middle of the layer respectively. Once the 
dimensionless numbers together with other physical input values are determined, the profiles of p a , 
T a , Pa, and S a can be evaluated. 


It follows from (23) that setting n = 1/(7 — 1) yields an adiabatic stratification ( dS a /dr = 0). 


We do so here so that the a subscript denotes a polytropic, adiabatic stratification. This provides an 
excellent first approximation to a stellar or planetary convection zone which is very nearly adiabatic 
due to the high efficiency of the convection. Though they cannot be strictly adiabatic because they 
must satisfy the Schwarzschild criterion ( dS/dr < 0), they are nearly adiabatic in the sense that 


d ,dS. 


(28) 


The bar denotes an average over latitude and longitude. Equation (28) is the basis of the so-called 


anelastic approximation in which the equations of motion are derived as perturbations about a 
static, often (but not necessarily) adiabatic, reference state f3T, 30] [32] . 

Although we do not employ the anelastic approximation here, an adiabatic, polytropic reference 
provides a useful starting point for setting up the initial conditions. However, the process cannot 


end here because this polytropic solution does not in general satisfy Eq.(21) and, since it does not 
satisfy the Schwarzschild criterion, it will not excite convection. 

4-3. Almost Flux Balance Approach 

In anelastic systems, the superadiabatic component of the entropy gradient ( dS/dr < 0) is 
assumed to be small so it can be specified indepedently of the adiabatic reference state [TO]. This 
amounts to setting po = p a , Tq = T ai and Pq = P a , and then solving equation (21) for dSo/dr : 


dSo 

dr 




1 , L , r, dT a ^ 

n T n / 47 rr 2 + PaGp ‘ r dr 


(29) 


where L is the luminosity. 

This procedure breaks down in fully compressible systems because equation (23) will only be 
satisfied to lowest order in e. This is a high price to pay merely to satisfy equation (21) which 


should have little bearing on the final dynamical equilibrium achieved after the onset of convection. 


It is more essential to satisfy the constitutive equations (22) and (231 precisely, together with the 


hydrostatic balance equation ( 20 ) to avoid a rapid initial restratification. 


We achieve this by introducing an extra step in the initialization procedure. As in anelastic 
systems, we compute the polytropic, adiabatic stratification as described in section |4.2| and we 
calculate T as defined in Eq.(29). However, unlike anelastic systems, we treat T as a target entropy 


gradient and then solve equations (20), (22), and (23) precisely using a separate finite difference 
code. In particular, we solve the following two equations for po and Pq using p a and P a as an initial 
guess 

r , 1 dp 0 gp 0 dp 0 

7 - = + and ~TZ = -Po9- 


c„ 


>0 dr 7Po' 


dr 


(30) 
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The temperature To is then given by (22). This process produces a superadiabatic, hydrostatic, 
splrerically-symmetric initial state that satisfies equations (20), (22), and (23). The entropy gradient 
will be equal to T but the thermal energy balance equation (21) will only be satisfied to lowest order 
in e. 

An alternative approach would be to solve all four equations (20)-(231 simultaneously for the 
four unknowns Pq, poj To, and Sq. Our Almost Flux Balance approach is much easier to implement 
and provides an effective way to initiate convection. 


5. Boundary Conditions 


The inner and outer boundaries are assumed to be impenetrable and free of viscous stresses 


^ = #(—> = t-(-) = °- 


(31) 


where V r , Vg, are the velocity components in spherical coordinates (r, 9, <f>). 

In addition, a constant heat flux f • r = L/(4irRf) is imposed at the bottom boundary and the 
temperature is fixed at the top boundary. 

The present CHORUS code employs 20 nodes including 8 corner points and 12 mid-edge points 
for each hexaliedral element in the physical domain. A precise treatment of the curved top and 
bottom boundaries of the spherical shells is then assured by the iso-parametric mapping procedure 
mentioned in Section 3. For calculating inviscid fluxes on element interfaces, an approximate 
Riemann solver |24j is generally used. However, exact inviscid fluxes are employed on top and 
bottom boundaries by using the fact that V r is precisely zero on spherical shell boundaries. A careful 
transformation between Cartesian and Spherical coordinate systems is conducted for computing 
viscous fluxes on two boundaries in order to ensure the stress-free conditions. 


6. Angular Momentum Conservation 

The equations of motion (l)-(3) express the conservation of mass, energy, and linear momentum. 
The conservation of angular momentum follows from these equations and the impenetrable, stress- 
free boundary conditions discussed in section [5j These are hyperbolic equations and we express 
them in conservative form when implementing the numerical algorithm, as disscussed in section [3] 
This, together with the spectral accuracy within the elements and the approximate Riemann solver 
employed at cell edges, ensures that the mass, energy, and linear momentum are well conserved as 
the simulation proceeds. However, we do not explicitly solve a conservation equation for angular 
momentum. Numerical errors including both truncation and round-off errors can result in a small 
change in angular momentum over each time step. Though these changes may be small, even a 
highly accurate algorithm can accummulate errors over thousands or millions of time steps that can 
compromise the validity of a simulation. This can be an issue in particular for unstructured grid 
codes like ours that solve the conservation equations in Cartesian geometries that are mapped to 
conform to the spherical boundaries. However, conservation of angular momentum can be violated 
even in highly accurate pseudo-spectral simulations, as reported by Jones et al. m- For this reason, 
we introduce an angular momentum correction scheme similar to one of the schemes described by 

Jones et al D3H- 

Two correction steps are taken in the CHORUS code to maintain constant angular momentum 
over long simulation intervals including: 
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Step 1. Calculate three Cartesian components of angular momentum explicitly, namely ( L x , L y , L z ). 
Step 2. Add a commensurate rigid body with rotating rate (5ft x , 5ft y , 5ft z ) to remove the angular 
momentum discrepancy. 

In Step 1, three Cartesian components are evaluated over the full spherical shell V as 

L x = pr(w sin 9 sin </> — v cos 9)dV, 

Jv 

L y = pr(u cos 9 — w svc 9 cos (j>)dV, (32) 

Jv 

L z = pr(v sin 9 cos <p — u sin 9 sin <p)dV. 

Jv 


Note that all three components are initially zero relative to the rotating coordinate system. 
The introduced rigid body rotating rate (5ft x ,5ft y ,5ft z ) in Step 2 is determined by 

5ft x = l x /i x , sn y = L y /i y , 5ft z = l z /i z , 


(33) 


where I x ,I y ,I z are the moment of inertia of the the spherical shell in the x , y and z directions 
respectively, and I x = f v p(y 2 + z 2 )dV , I y = f v p(x 2 + z 2 )dV, and I z = f v p(x 2 + y 2 )dV. Once 
(5ft x , 5Cl y , Stt z ) is obtained, the Cartesian velocity components u, v, and w will be updated at each 
solution point using 


(34) 


Unew = u 0 id + 5ft, z \Jx 2 + y 2 sin(atan2(j/, x)) — 5ft y \Jx 2 + z 2 cos(atan2(x, z)), 

Vnew = Void ~ 5ft z \/x 2 + y 2 cos(atan2(?/, x)) + 5ft x ^/y 2 + z 2 sin(atan2 (z,y)), 
uinew = Wold + 5ft y sJx 2 + z 2 sin(atan2(x, z)) — 5ft x yfy 2 + z 2 cos(atan2(z, y)), 

where function atan2(mi,m 2 ) = 2arctan( , mi = -). 

We note that this correction procedure applies equally well in the case of an oblate, rapidly- 
rotating star. The total moment of inertia in each direction is numerically calculated by summing 
up the moment of inertia in each element of the unstructured grid so it works regardless of the 
shape of the star. 

This correction procedure is computationally expensive if performed at every time step. Thus, 
in all CHORUS simulations, we only do correction every 5000 time steps. Fig|4] shows the time 
evolution of 6ft X) 5ft y and 5ft z from a representative CHORUS simulation. This demonstrates 
the high accuracy of the numerical algorithm since the relative angular momentum error never 
exceeds 2 x 10~ 6 even at this 5000-step correction interval. Furthermore, it demonstrates that the 
cumulative long-term errors in the angular momentum components are well controlled. 


7. Code Verification 

We verify the CHORUS code by comparing its results to the well-established Anelastic Spherical 
Harmonic (ASH) code I.'SIT 1'?. .17 . However, we acknowledge that this comparison is not ideal since 


the two codes solve different equations. As mentioned in section 4.2 the equations of motion in the 


anelastic approximation are obtained by linearizing the fully compressible equations (l)-(3) about 
a hydrostatic, spherically-symmetric reference state. So we would only expect CHORUS and ASH 
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Figure 4: D eviations from angular momentum conservation expressed in terms of angular velocity variations according 
to Eq. ( |.13[ i. 6 fl x , Sfly, and Sfl z are sampled every 5000 iterations just before the correction procedure and normalized 
to the rotation rate of the coordinate system Hq. 


to agree in the limit e —> 0, where e is the normalized radial entropy gradient defined in Eq. (28). 
A thorough comparison between the two systems would involve a linear and nonlinear analysis 
demonstrating convergence as e —» 0. This lies outside the scope of the present paper. Here we 
focus on defining two benchmark simulations, patterned after the gaseous atmosphere of Jupiter 
and the convective envelope of the Sun, and compare the results from CHORUS and ASH. Despite 
subtle differences in the model equations and substantial differences in the numerical method, we 
demonstrate good agreement between the two codes. This serves to verify CHORUS and to pave 
the way for future applications. 


7.1. The ASH Code and Anelastic Benchmarks 

ASH is a multi-purpose code designed to simulate the hydrodynamics (HD) and magnetohydro¬ 
dynamics (MHD) of solar and stellar interiors in global spherical geometries. It was first developed 
over 15 years ago and has remained at the leading edge of the field ever since, continually improving 
in its physical sophistication and parallel efficiency on high-performance computing platforms. ASH 
results have appeared in over 100 publications with applications ranging from convection and dy¬ 
namo action in solar-like stars and fully-convective low-mass stars, to core convection and dynamo 
action in massive stars, to MHD instabilities and stably-stratified turbulence, to the generation of 
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and transport by internal gravity waves, to tachocline confinement, to flux emergence, to the HD 
and MHD of red giants 01331351 EZ| . 

ASH is based on the anelastic approximation (sec. 4.2) and uses a poloidal-toroidal decomposi¬ 
tion of the mass flux to ensure that the anelastic form of the mass continuity equation is satisfied 
identically (V-(po u ) =0). It is a pseudo-spectral code that uses triangularly-truncated spherical 
harmonic basis functions in the horizontal dimensions. Although earlier versions of ASH employed 
Chebyshev basis functions in the radial dimension, the version presented here uses a centered, 
fourth-order finite difference scheme that was introduced to improve parallel efficiency. The radial 
grid is uniformly spaced for the simulations presented here and the boundary conditions are as 
specified in section [5] Time stepping is accomplished using an explicit Adams-Baslrforth scheme 
for the nonlinear terms and a semi-implicit Crank-Nicolson scheme for linear terms, both second 
order. 

ASH is one of four global anelastic codes that were validated using a series of three carefully 
defined benchmark simulations presented by Jones et al.pJ)]. All three benchmarks had shell ge¬ 
ometries and dimensional parameters that were chosen to represent deep convection in Jupiter’s 
extended atmosphere but they differed in their degree of magnetism (one was non-magnetic) and 
turbulent intensity (spanning laminar and turbulent dynamo solutions). In order to concisely rep¬ 
resent the effective parameter space, the benchmark simulations were specified through a series of 
non-dimensional parameters as defined in Eq.(24), namely R a , P r , Ek , and /3. These benchmarks 
made use of a hydrostatic, adiabatic, polytropic reference state as described in sec. |4.2[ specified 


through the additional nondimensional parameters N p and n (the polytropic index). All four codes 
used similar numerical methods (pseudospectral, spherical harmonic) and agreed to within a few 
percent for a variety of different metrics of physical quantities. 

The benchmarks we define here are inspired by the anelastic convection-driven dynamo bech- 
marks of Jones et al. [lO]. However, some modifications to the anelastic benchmarks are necessary in 
order to ensure that they are consistent with the fully compressible equations solved by CHORUS. 
We already discussed one example of this when defining the initial conditions in sec. |4.3| 

Another significant modification of the anelastic benchmarks that we introduce here concerns 
the Mach number. An implicit requirement of the anelastic approximation is that the Mach number 
of the flow is much less than unity M a = U/c s « 1, where c s is the sound speed. This is well 
justified in stellar convection zones where M a is typically less than 1CU 3 . In a fully compressible 
code this places a severe constraint on the allowable time step permitted by the Courant-Freidrichs- 
Lewy (CFL) condition At < S/c s where S is some measure of the minimum grid spacing. Anelastic 
codes are not subject to this constraint. If the Mach number is low, the CFL constraint imposed 
by the sound speed is much more stringent than that imposed by the flow field At < S/U. In the 
future we will mitigate this constraint through the use of implicit and local time stepping. Here 
we address it by defining benchmark problems for which the Mach number is low (justifying the 
anelastic approximation in ASH) but not too low (mitigating the CFL constraints of CHORUS). 

The CFL constraint arising from the sound speed can be a major issue for global, rotating 
convection simulations where the equilibration time scale is much longer than the dynamical time 
scale. If one neglects structural stellar evolution, the longest time scale in the system is the thermal 
relaxation timescale T r = E /X, which can exceed 10 5 years in stars; by comparison the dynamical 
time scale is of order one month. However, a more relevant time scale for equilibration of the 
convection is the thermal diffusion time scale Td = cP/k, which is ~ 65.74 days for the Jupiter 
benchmark and ss 58.44 days for the solar benchmark (T r ~ 3.55 x 10 4 days and 5.45 x 10 4 days 
respectively). 


15 




7.2. Metrics of Physical Quantities 

Before proceeding to the simulation results, we first define several important metrics that provide 
a means to compare CHORUS and ASH. These include the mean kinetic energy relative to the 
rotating reference frame and its mean-flow components, namely the the differential rotation (DRKE) 
and the meridional circulation (MCKE): 


KE=yJ v lp(u.n)dV, 

DRKE = f ^p(U^,) 2 7' 2 sin ddrdddcj), 
V Jv 2 


MCKE = 


b 


i 


p({V r } 2 + ( Vg) 2 )r 2 sin ddrdddcj) 


(35) 

(36) 

(37) 


v 


where V denotes the volume of the computational domain and angular brackets denote averages 
over longitude. 

The growth rate of the kinetic energy is defined as 

d(lnKE) 


a = 


dt 


(38) 


From Eq. four components of the energy flux are involved in transporting energy in the radial 
direction, namely the enthalpy flux F e . kinetic energy flux Fk, radiative flux F r and entropy flux 
F u . In a statistically steady state, these four fluxes together must account for the full luminosity 
imposed at the bottom boundary: 


where 


+ Fk + F r + F u — -F* — ^ r 2 > 

(39) 

F e = pC p V r (T — T), 

(40) 

Fk = -pV r (u-u), 

(41) 

dT 

F r = K r pCp-^—^ 

(42) 

P -ypdS 

F u = ~k P T—. 

(43) 


and 


p , T, and S denote the mean density, temperature and entropy, averaged over horizontal surfaces. 
On each horizontal surface, the mean Mach number is defined as 


where 


V„ 


and the mean sound speed is 


Ma = 


V„ 


c, 



(44) 


(45) 

(46) 
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Dimensionless parameters 

Ek = 10 -3 , N p = 5 , /3 = 0.35, Ra = 351,806, Pr = 1, n = 2.0 
Defining physical input values 

R 0 = 7 x 10 9 cm, n 0 = 1.76 x HT 4 s" 1 , M = 1.9 x 10 3 ° g, Pi = 1.1 g cm~ 3 
m = 3.503 x 10 7 erg g~ l K~\ G = 6.67 x lO" 8 g~ l cm 3 s" 2 , n r = 0 

Derived physical input values 

Ri = 2.45 x 10 9 cm, d = 4.55 x 10 9 cm, v = 3.64364 x 10 12 cm 2 s" 1 , 
k = 3.64364 x 10 12 cm 2 s _1 , 7 = 1.5 

Other thermodynamic quantities 

L = 7.014464 x 10 32 erg s" 1 , C p = 1.0509 x 10 s erg g~ 1 K~ x 


Table 1: Parameters for the Jupiter benchmark 


7.3. Jupiter Benchmark 

The deep, extended outer atmosphere of Jupiter is thought to be convectively unstable [38] . 
This has motivated substantial work on the internal dynamics of giant planets and inspired the 
parameter regimes chosen for the anelastic benchmarks of Jones et al. m- Our first benchmark 
for comparing CHORUS and ASH is similar to the hydrodynamic benchmark of Jones et al. uni 
apart from the thermal boundary conditions. Whereas we impose a fixed heat flux on the lower 
boundary and a fixed temperature on the upper (sec. [5]), Jones et al. fix the specific entropy S on 
both boundaries. Note that this means that the Rayleigh number R a , defined in terms of AS in Eq. 
(24), does change somewhat in our simulations as the convection modifies the entropy stratification. 
This is in contrast to Jones et al. where it was held fixed. 

The parameters for this case are specified in Table [l] The value of R a listed is the initial value, 
before convection ensues. The number of DOFs used for the CHORUS simulation is about a factor 
of five larger than that used for the ASH simulation, as indicated in Table [3] Both simulations use 
the same boundary conditions and initial conditions, apart from the random perturbations needed 
to excite convection which are generated independently by each code. 


7.3.1. Exponential Growth and Nonlinear Saturation 

The evolution of the kinetic energy densities for both CHORUS and ASH simulations are illus¬ 
trated in Fig |5(a)| As mentioned above, both simulations start with the same background stratifi¬ 
cation but with different random perturbations. The small amplitude of the initial perturbations 
ensures that each simulation begins in the linear regime. For each simulation there is an initial 
adjustment period before the flow locks on to the fastest-growing eigenmode which then grows 
exponentially. 

The initial establishment period is different for each simulation, lasting roughly 8 days for the 
CHORUS simulation and 2 days for the ASH simulation. However, this is to be expected from the 
different mix of random perturbations. A meaningful comparison can only be made between the 
two codes after they reach the exponential growth phase, which is followed by nonlinear saturation 
and a subsequent equilibrium. In other words, the two features of Fig. |5(a)~| that should be compared 
are the slope of the line in the linear growth phase and the value of the kinetic energy density after 
each simulation has saturated and equilibrated. 
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(a) 


(b) 


Figure 5: (a) Kinetic energy density and (6) growth rate for the Jupiter benchmark. 


This first point of comparison, namely the growth rate, is highlighted in Fig |5(b)| The peak 
values achieved by each simulation in the linear regime reflect the growth rate of the preferred linear 
eigenmode and agree to within about 2%; cr = 2.38 x 10 -5 for CHORUS and a = 2.33 x 10 -5 for 
ASH. We define the nonlinear saturation time T n as the time at which the growth rate first crosses 
zero after the exponential growth phase. 

Again, the saturation time is different between the two simualtions because of the random nature 
of the initial conditions. Thus, in order to compare the two cases in the nonlinear regime we define 
a sampling time T s to be 10 days after the saturation time, T n . This is T s = 23.51 days for the 
CHORUS simulation and T s = 21.15 days for the ASH simulation. Averaged between (T s — 2) days 
and T s days, KE = 4.33 x 10 5 erg cm~ 3 for the CHORUS simulation and I\E = 4.15 x 10 5 erg 
cm~ 3 for the ASH simulation. The difference is about 4%. Another point of comparison is the 


relative magnitude of the mean flows, as quantified by the DRKE and MCKE defined in sec. 7.2 
At the sampling time T s , DRKE/KE = 0.113 and MCKE/KE = 2.45 x 10" 4 for the CHORUS 
simulation while DRKE/KE = 0.134 and MCKE/KE = 2.58 x 10 -4 for the ASH simulation. This 
corresponds to a difference of about 16% and 5% for the DR and MC respectively. Small differences 
between the two codes of order several percent are to be expected due to differences between the 
compressible and anelastic equations. Furthermore, the relatively large difference in the DRKE/KE 
likely comes about because the simulations are not fully equilibrated. We address these issue further 


in the following section (sec. 7.3.21. 


7.3.2. Mach Number, e and Flux Balance 

As mentioned in secs. |4.2| and |7.1[ the comparisons between ASH and CHORUS are only mean¬ 
ingful if the stratification is nearly adiabatic (e << 1) and the Mach number is small. These 
conditions are required for the validity of the anelastic approximation. Fig. |6(a)| demonstrates that 
these conditions are met, but that departures are significant. In particular, we expect that the 
anelastic and compressible equations are only equivalent to lowest order in e, which reaches a value 


as high as 0.07 in the upper convection zone (Fig 6(b)). 

In Fig|6( a)j the mean Mach number is minimum at the bottom and increases with radius and 
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reaches the maximum (about 0.012) at the top in both the CHORUS and ASH simulations. As 
defined in Eq.(28), e is proportional to the mean entropy gradient (pG Identical initial entropy 
gradients for the CHORUS and ASH simulations implies that they have the same initial degrees 
of adiabaticity. When convection is present, the associated energy flux leads to a redistribution of 
entropy, tending to smooth out the entropy gradient. This is the reason that e in Fig |6(b)| becomes 
smaller than the initial near r = 0.92 R 0 . 

The location r = 0.92i?. o corresponds to where the efficiency of the convection peaks. This is 
demonstrated in Fig[7] which shows the components of the energy flux defined in sec. |7.2[ The 
almost flux balance initialization described in sec. |4.3| establishes an entropy stratification that 


carries most of the energy flux through entropy diffusion F u « F*, with a slight over-luminosity of 
about 7% in the upper convection zone. Once convection is established, it carries roughly 15% of 
this flux, flattening the entropy gradient and reducing F u . The kinetic energy flux F\ for this case 
is negligible. The maximum values of F e at r = 0.92I? o for the CHORUS and ASH simulations are 
0.1536 and 0.1484 respectively, showing the difference of 3.5%. 




r/R 0 r/R 0 


(a) 


(b) 


Figure 6: (a) Mach number and (b) e for the Jupiter benchmark, plotted at the sampling time T s . In (b), the initial 
condition is also plotted for comparison (dash-dotted line). The CHORUS and ASH curves are nearly indistinguish¬ 
able. 


In a nonlinear equilibrium state, the sum of the normalized fluxes in Fig.[7]should be unity. This 
is clearly not the case; as mentioned above both simulations are over-luminous by about 7% due 
to the initial entropy stratification. This will eventually subside but the process is slow, occurring 
gradually over a time scale that is longer than the thermal diffusion time scale of Td ~ 65.7 4 days 
but shorter than the thermal relaxation time scale of T r ~ 3.55 x 10 4 days (see sec. 7.1). This 
is demonstrated in Fig(8] which shows the flux balance in the ASH simulation after 1800 days. 
Given the greater computational cost of CHORUS (sec. 7.5), and the satisfactory agreement at the 
sampling time (within the expected order e), we choose not to run the CHORUS simulation to full 
equilibration. 


7.3.3. Convection Structure 

The structure of the convection at the sampling time is illustrated in Fig. [9] Though this is 
well into the nonlinear regime, both simulations are dominated by a series of columnar convective 
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Figure 7: Components of the radial energy flux for the Jupiter benchmark integrated over horizontal surfaces at 
the sampling time T s for (a) CHORUS and ( b ) ASH. All values are normalized by the total flux F * imposed at the 
bottom boundary. 



Figure 8: Components of the radial energy flux as in Fig. [7] but at a much later time, t ~ 1800 days, computed with 
the ASH code, showing an equilibrium solution. 
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rolls approximately aligned with the rotation axis but sheared slightly in the prograde direction at 
low latitudes by the differential rotation. These are the well-known ‘banana cells’ characteristics 
of convection in rotating spherical shells and most apparent for laminar parameter regimes [5]. 
Though this is well into the nonlinear regime, they reflect the preferred linear eigenmodes and are 
well described by a single sectoral spherical harmonic mode with £ = m, where £ and to are the 
spherical harmonic degree and order. The degree and order £ and to can also be interpreted as the 
total wavenumber and the longitudinal wavenumber respectively. 

Close scrutiny of Fig. [9] reveals that the two simulations exhibit a slightly different mode struc¬ 
ture, with CHORUS selecting an to = 20 mode and ASH selecting an m = 19 mode. This is 
demonstrated more quantitatively in Fig jTo] which shows the spherical harmonic spectra of the 
velocity field on the horizontal surface r = 0.98i? o as a function of spherical harmonic degree £ 
(summed over to) at the sampling time T s . In the CHORUS simulation, the spectra of radial ve¬ 
locity V r , meridional velocity Vg. and zonal velocity peak at £ = 20 and its higher harmonics, 
£ = 40, 60 and 80. In comparison, the velocity spectra in the ASH simulation peak at £ = 19, 38, 
57 and 76. 

This level of agreement is consistent with the expected accuracy of the anelastic and compressible 
systems. As discussed by Jones et al. uni, the linear growth rates of the to = 19 and to = 20 modes 
in this hydrodynamic benchmark are the same to order e and even a single anelastic code may 
choose one or the other depending on the random initial perturbations. The laminar nature of 
this benchmark highlights these small differences; in more turbulent parameter regimes where the 
Rayleigh number far exceeds the critical value, a broad spectrum of modes is excited and the 
results are less sensitive to the details of the initial conditions and the linear eigenmodes. This is 
demonstrated by the solar benchmark described in sec. |7.4| which is in a more turbulent parameter 
regime and which exhibits closer agreement between the velocity spectra. 


7.3-4- Mean Flows 

The mean (averaged) flows for the Jupiter benchmark are shown in Fig. 11 All averages span 
2 days (about 4.8 rotation periods), starting at the sampling time T s . 

The meridional circulation is expressed in terms of a stream function, T, defined as 

_ 194/ _ 94/ 

r sin 9(pV r ) = and r sin 9 (pVg) = — , 

r 06 or 

and the differential rotation is expressed in terms of the angular velocity 

n = ^-(^- + n 0 ). 

Z7r r sin 0 


(47) 


(48) 


We define thermal variations S and T by averaging over longitude and time and then subtracting 
the spherically-symmetric component [£ = m = 0) in order to highlight variations relative to the 
mean stratification. 

The CHORUS and ASH results in Fig. 0 correspond closely with a few notable exceptions. 
Near the equator in the upper convection zone, SI in Fig|ll[d) is somewhat smaller than that 


FigJTTJ A,). This is also reflected by the lower DRKE/KE noted in sec. 7.3.1 As mentioned 


there, this discrepancy may in part be because the simulations are not strictly in equilibrium at 
the sampling time (Fig. [7]). We would expect the correspondence to improve if we were to run 
CHORUS for several thousand days, giving the mean flows ample time to equilibrate along with 
the stratification. The wiggles in the S' plot for CHORUS (Fig. [TT|)fc)) can be attributed to two 
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(a) CHORUS simulation 
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(b) ASH simulation 


Figure 9: Mollweide projections of radial velocity V r for the Jupiter benchmark at the horizontal surface r = 0.95R o , 
taken for (a) CHORUS and (6) ASH at the sampling time T s . Red and blue tones denote the upflow and downflow 
as indicated by the color bar. 
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(a) 


(b) 


(c) 


Figure 10: Shown are the power spectra of (a) radial velocity V r , (6) meridional velocity V$, and (c) zonal velocity 
as a function of spherical harmonic degree l for the Jupiter benchmark at r = 0.98_R o and t = T s . 


factors. First, unlike ASH, the relevant state variable in CHORUS is the total energy E. The specific 
entropy must be obtained from E by subtracting out contributions from the kinetic energy and the 
mean stratification to obtain the pressure variations, from which S is computed (using also the 
density variation). This residual nature of S 1 along with the post-processing step of interpolating 
the CHORUS results onto a structured, spherical grid both contribute numerical errors. In ASH, 
by contrast, S' is a state variable. The correspondence of Figs. &) and (/) despite these numerical 
errors is a testament to the accuracy of CHORUS. 

7-4- Solar Benchmark 

Having gained confidence in simulating the Jupiter benchmark, we now look into defining a 
benchmark that has more in common with the Sun. This serves two purposes. First, it helps 
verify the CHORUS code by probing a different region in parameter space, this one more turbulent. 
Second, it makes explicit contact with solar and stellar convection which is the primary application 
we had in mind when developing CHORUS. 

The most significant differences between the solar and the Jupiter benchmarks include the 
Rayleigh number R a , the density stratification N p , and the radiative heat flux F r . The value of 
R a is about four times larger in the solar benchmark. This, combined with the larger value of E^ 
implies that the flow is more turbulent and less rotationally constrained (larger Rossby number). 
In the middle of the convection zone, the Rossby number ~ 1.56 x 1CU 2 for the Sun benchmark and 
it is larger than the Jupiter benchmark ~ 3.81 x 10 -4 . As in the Sun, the radiative flux F r (which 
was set to zero in the Jupiter benchmark) carries energy through the bottom boundary, extending 
into the lower convection zone (see Fig.[l4|. This allows us to set the radial entropy gradient dS/dr 
to zero at the lower boundary, which is what marks the base of the convection zone in the Sun. 
Another significant difference is the aspect ratio /3, for which we use a solar-like value of 0.73. 

As mentioned in secs. [l]and |7.1| a challenge in simulating solar convection with a compressible 
code like CHORUS is the low value of the Mach number. We address this challenge by scaling up 
the luminosity L by a factor of 1000 relative to the actual luminosity of the Sun. According to 
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Figure 11: Mean flows in the Jupiter benchmark from CHORUS (top row) and ASH (bottom row). Shown are (a, 
e) meridional circulation expressed in terms of the stream function with red and blue denoting clockwise and 
counterclockwise circulations respectively, (6, /) specific entropy perturbation S , (c, g) temperature perturbation 
T , and ( d , h) differential rotation, expressed in terms of the angular velocity fl. The top and bottom plots in each 
column share the same color bar. 


24 











































Dimensionless parameters 

E k = 2.447 x 10” 3 , N p = 3, /3 = 0.736762, R a = 1,428,567, Pr = 1, n = 1.5 
Defining physical input values 

R 0 = 6.61 x 10 10 cm, D 0 = 8.1 x 10” 5 s” 1 , M = 1.98 891 x 10 33 g, = 0.21 g cm~ 3 
91 = 1.4 x 10 s erg g~ x K~\ G = 6.67 x 10” 8 g~ l cm 3 s~ 2 

Derived physical input values 

Ri = 4.87 x lO 10 cm, d = 1.74 x 10 10 cm, v=6.0 x 10 13 cm 2 s” 1 , 
k = 6.0 x 10 13 cm 2 s” 1 , 7 = 5/3 

Other thermodynamic quantities 

L = 3.846 x 10 36 erg s”\ C p = 3.5 x 10 8 erg g~ x K~ l 


Table 2: Parameters for the solar benchmark 


the mixing-length theory of convection, the rms velocity should scale as L 1 / 3 . Thus, to preserve a 
solar-like value of the Rossby number (important for achieving reasonable mean flows), this suggests 
that we would need to scale up the rotation rate f2 0 by a factor of approximately 10 relative to the 
actual Sun. We then chose values of of k and v to give practical values for R a and E k . Table [2] 
gives the parameters for the solar benchmark. 

The radiative flux F r is parameterized by expressing the radiative diffusion as n r = A(co + cio;-|- 
c 2 w 2 ), where c 0 = 1.5600975 x 10 s , ci = -4.5631718 x 10 7 , c 2 = 3.3370368 x 10 6 , and w = r x 10” 10 . 
The parameter A is chosen so that F r = L/(4nr 2 ) on the bottom boundary. As for the Jupiter 
benchmark, we define a sampling time T s that is 10 days after the nonlinear saturation time. The 
sampling times are T s = 15.10 days and T s = 14.59 days for CHORUS and ASH respectively. The 
number of DOFs used for the CHORUS simulation is about a factor of 6.8 larger than that used 
for the ASH simulation, as indicated in Table [3] 


7.4-1. Exponential Growth and Nonlinear Saturation 

The volume-averaged kinetic energy densities for both the CHORUS and ASH simulations are 
plotted in Fig[i~2|(a). As in the Jupiter benchmark, the random initial conditions lead to differences 
in the initial transients but once the preferred linear eigenmode begins to grow the two simulations 
exhibit similar growth rates and nonlinear saturation levels. Averaged between (T s —1) and T s days, 
KE = 9.02 x 10 8 erg cnrT 3 for the CHORUS simulation and KE = 8.58 x 10 8 erg cm ~ 3 for the 
ASH simulation. The difference is about 5.13%. In the linear regime, they grow very fast and the 
ranges of a well-defined linear unstable regime are narrow as shown in Fig[l2|(b) for both CHORUS 
and ASH simulations. Thus, only their maximum growth rates in that regime are compared. A 
difference of only 4% is present between a max = 6.02 x 10 -5 of the CHORUS simulation and 
o'max = 6.29 x 10~ 5 of the ASH simulation. At the sampling time T s , DRKE/KE= 4.66 x 10” 1 
and MCKE/KE= 9.02 x 10” 4 for the CHORUS simulation while DRKE/KE= 5.08 x 10” 1 and 
MCKE/KE= 1.07 x 10” 3 for the ASH simulation (a discrepancy of about 9% and 19% respectively). 
As in the Jupiter benchmark, this relatively large difference in mean flows is likely attributed to the 
immature state of the simulations, which has not yet achieved flux balance by the sampling time 
(sec. 7.4.2). 
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Figure 12: (a) Kinetic energy density and (6) growth rate for the solar benchmark. 


7-4-2. Mach Number, e and Flux Balance 

As in the Jupiter benchmark, a low Mach number is also achieved in this solar benchmark, with 
the maximum Mach number around 0.017 and 0.015 for the CHORUS simulation and the ASH 
simulation, respectively, as shown in Fig 13(a) The value of e also peaks in the upper convection 
zone (Fig. |13(bjj ). Reflecting the flatter entropy gradient at the sampling time T s , e becomes smaller 
where the convective efficiency is largest, near r = 0.95i? o . 

The CHORUS and ASH simulations also have similar flux balances as shown in Fig fid} Both 
are over-luminous (F e + F/. + F r + F u > A*) due mainly to the large diffusive and entropy fluxes 
F u and F e , which have not yet equilibrated. This imbalance subsides by i ~ 540 days as verified by 
an extended ASH simulation. However, as with the Jupiter benchmark, we have not run CHORUS 
to full equilibration because of the computational expense. Increasing the luminosity further will 
mitigate this relaxation time and we plan to exploit this for future production runs. 

The radiative flux F r carries energy through the bottom boundary and dominates the heat 
transport in the lower convection zone while the entropy flux F u carries energy through the top 
boundary and the upper convection zone. The enthalpy flux F e gradually increases towards the top 
until peaks near r = 0.95i? o , and then drops down to zero rapidly as it approaches the impenetrable 
top boundary. At r = 0.95i? o , F e = 0.50 and F u = 0.78 from the CHORUS simulation while 
F e = 0.42 and F u = 0.76 from the ASH simulation. 


7-4-3. Convection Structure 

The structure of the convection in the solar benchmark is illustrated in FiglU One would 
not expect the instantaneous flow field in a highly nonlinear simulation to correspond exactly 
between two independent realizations but the qualitative agreement is promising. This qualitative 
agreement is confirmed quantitatively by comparing the velocity spectra in Fig |16| The radial 
velocity spectrum peaks at £ = 15 — 25 (Figfl6](a)) and the meridional velocity spectrum exhibits 
substantial power in two wavenumber bands, namely £ = 15 — 25 and £ = 35 — 45 (Figfl6|(fc)). For 
the zonal velocity spectra (FigJT6]( c)), three modes £ = 2,4, and 6 are prominent in the low £ (< 10) 
range. After peaking at £ = 15 — 25, the high-^ power decreases exponentially in amplitude. Some 
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Figure 13: (a ) Mac h number and ( b ) normalized entropy gradient e for the solar benchmark at the sampling time 
T s . As in Fig. |6(bJ[ the initial profile is included in frame (6), where the CHORUS and ASH curves at T s are slightly 
steeper and indistinguishable from one another. 



Figure 14: Components of the normalized, horizontally-integrated radial energy flux as in Fig. [7] but for the solar 
benchmark. Each snapshot corresponds to the sampling time T s . 
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of the small discrepancies between the two curves in each plot can likely be attributed to random 
temporal variations that would cancel out with some temporal averaging. 
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Figure 15: Mollweide projections of the radial velocity V r in the solar benchmark at r = 0.95 R 0 and t = T s for (a) 
CHORUS and ( b ) ASH. Red and blue tones denote the upflow and downflow as indicated by the color bar. 


7.4-4’ Mean Flows 

The meridional circulation from the CHORUS simulation (see Figfl7](a)) and the ASH simulation 
(see Fig© e)) exhibit similar flow patterns with most circulations concentrating between the low 
latitudes and middle latitudes, outside the so-called tangent cylinder, namely the cylindrical surface 
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(a) (b) (c) 


Figure 16: Power spectra of (a) radial velocity V r , (b) meridional velocity V#, and (c) zonal velocity Vq as in Fig. |10| 
but for the solar benchmark at r = 0.98 R 0 and t = T s . 


aligned with the rotation axis an d ta ngent to the base of the convection zone. T he s pecific entropy 


perturbation S are shown in Fig 17 b) for the CHORUS simulation and in Fig 17 f) for the ASH 


simulation. In both sim ulati ons, the contours of S are symmetric a bou t the equator. By comparing 


the contours of T in Fig 17 c) for the CHORUS simulation and Fig 17 g) for the ASH simulation, a 


good agreement is achieved. For both simulations, they also have similar differential rotation profiles 
as shown in Figfl7|)d) and FigJT7th). Some of the small discrepancies can likely be attributed to the 
flux imbalances shown in Fig |14| causing mean flows to vary slowly as the simulations equilibrate 
from different random initial conditions and nonlinear saturation states. Residual random temporal 
fluctuations may also be present despite the (short) time average, particularly for the meridional 
circulation which is a relatively weak flow with large fluctuations [5]. 


7.5. Code Performance 

As demonstrated in sec. [3j the CHORUS code achieves excellent scalability out to 12k cores. This 
is strong scaling for intermediate-resolution simulations. We expect higher-resolution simulations to 
scale even better to tens of thousands of cores. In this section we consider CHORUS’s performance 
relative to the anelastic code ASH. 

The computational efforts are summarized in Table [3] The resolution in ASH is expressed as 
N r x Ng x N ^ where N r , Ng, and N $ are the number of grid points in the radial, latitudinal, and 
longitudinal directions respectively. However, due to pseudo-spectral de-aliasing and the symmtery 
of the spherical harmonics, the effective number of DOFs for ASH is N r (l max + 1 )( m axi where 
(-max = (2-/Vg — l)/3 is the maximum spherical harmonic mode. For both benchmarks, Ng = 256 
and (max — 170. 

The total number of core hours needed to run 10 days is much larger for CHORUS than it 
is for ASH; by more than three orders of magnitude for the Jupiter benchmark and by a factor 
of 75 for the solar benchmark. Much of this is due to the smaller time step required by the 
compressible scheme and the higher number of degrees of freedom used in running the CHORUS 
benchmarks. Furthermore, the ASH simulations were run on only 71 cores whereas the CHORUS 
runs typically employed several thousand.Thus, imperfect scaling is a factor, as is a difference in 
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Figure 17: Mean flows as in Fig(TT] but for the solar benchmark. Top and bottom rows correspond to CHORUS 
and ASH respectively and all quantities are averaged over longitude and over a two-day time interval beginning at 
T s (spanning 2.2 rotation periods), (a, e) mass flux stream function 'F with red and blue denoting clockwise and 
counterclockwise circulations respectively, (6, /) specific entropy perturbation S , (c, g) temperature perturbation 
T , and (d, h) angular velocity profile. The top and bottom plots in each column share the same color bar. 
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Jupiter Case 

Sun Case 

code 

CHORUS 

ASH 

CHORUS 

ASH 

Resolution 

294,912 elements 

129 x 256 x 512 

307,200 elements 

100 x 256 x 512 

DOFs 

18,874,368 

3,750,030 

19,660,800 

2,907,000 

Time step (s) 

1.5 

533 

4 

20 

Core hours per time step 

7.53 x 10 -2 

7.03 x 10 -3 

8.03 x 10 -2 

5.32 x 10“ 3 

Number of iteration 

required to run 10 days 

576,000 

1,621 

216,000 

43,200 

Number of core hours 

needed to run 10 days 

43,380 

11.4 

17,352 

230 


Table 3: Computational efforts for the Jupiter and solar benchmarks. 


the computational platform. The ASH simulations were run with the Intel Xeon E5-2680v2 (Ivy 
Bridge) cores on NASA’s Plieades machine (2.8 GHz clock speed, 3.2GB/core memory) whereas 
the CHORUS simulations were run with the Intel Xeon E5-2670 (Sandy Bridge) cores on NCAR’s 
Yellowstone machine (2.6 GHz clock speed, 2 GB/core memory). Furthermore, CHORUS uses 
a fourth-order accurate five-stage explicit Runge-Kutta method [55] whereas ASH uses a simpler 
second-order mixed Adams-Bashforth/Crank-Nicolson time stepping. This also contributes to the 
larger number of core hours per time step used by CHORUS (Table [3]). 

Though ASH out-performs CHORUS for these simple benchmark problems, it must be remem¬ 
bered that such problems are ideal for pseudo-spectral codes; relatively low resolution, laminar runs 
dominated by a limited number of spherical harmonic modes. The real potential of CHORUS will 
be realized for high-resolution, turbulent, multi-scale convection where its superior scalability and 
variable mesh refinement will prove invaluable. It can also be used for studying physical phenomena 
such as core convection and oblate stars that are challenging or even inaccessible to codes that use 
structured, spherical grids. Furthermore, there is much potential for improvement in the efficiency 
of CHORUS; we intend to implement an implicit time marching scheme and a p-multigrid method 
[40] as well as local time stepping in the future, and to optimize the numerical algorithm for higher 
performance on heterogeneous (CPU/GPU) architectures that are very suitable for data structures 
of the SDM. 

8. Summary 

We have developed a novel high-order spectral difference code, CHORUS, to simulate stellar 
and planetary convection in global spherical geometries. To our knowledge, the CHORUS code is 
the first stellar convection code that employs an unstructured grid, giving it unique potential to 
simulate challenging physical phenomena such as core convection in high and low-mass stars, oblate 
distortions of rapidly-rotating stars, and multi-scale, hierarchical convection in solar-like stars. 

The CHORUS code is fully compressible, which gives it advantages and disadvantages over 
codes that employ the anelastic approximation. On the one hand, the hyperbolic nature of the 
compressible equations promotes more efficient parallel scalability over the (elliptical) anelastic 
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equations. Indeed, we demonstrated that the CHORUS code does achieve excellent strong scala¬ 
bility for intermediate-size problems extending to 12,000 cores. We expect even better scalability 
for higher-resolution problems. Furthermore, the fully compressible equations are required to ac¬ 
curately capture the small-scale surface convection in solar-like and less massive stars where Mach 
numbers approach unity and where the anelastic approximation breaks down. On the other hand, 
the CFL constraint imposed by acoustic waves places strict limits on the allowable time step for 
simulating deep convection in most stars and planets, where the Mach number is much less than 
unity. We intend to address this constraint in the future by implementing implicit and local time 
stepping schemes. 

To verify the CHORUS code, we defined two benchmark simulations designed to bridge the gap 
between fully compressible and anelastic systems. This allowed us to compare CHORUS results 
to the well-established ASH code which employs the anelastic approximation. The two benchmark 
cases were formulated to simulate convection in the Jupiter and Sun. Metrics of physical quantities 
sensitive to the convective driving and structure such as the linear growth rate and the total KE 
agree to lowest order in e, the stratification parameter upon which the validity of the anelastic 
approximation is based. Mean flows exhibit larger variations (5-20%) that may be attributed to 
flux imbalances at the sampling time. Better agreement between the two codes could likely be 
achieved by running the simulations longer but full equilibration would require at least an order of 
magnitude more computing time. We did not believe this computational expense was warranted 
given the good agreement at the sampling time. The level of agreement is remarkable considering 
that the CHORUS and ASH codes not only solve different equations (compressible versus anelastic) 
but also they employ dramatically different numerical algorithms. Thus, we consider the CHORUS 
code verified. 

Future applications of CHORUS will focus on the applications discussed in §1, namely the inter¬ 
action of convection, differential rotation, and oblateness in rapidly-rotating stars, core convection 
in high and low-mass stars, hierarchical convection in solar-like stars, and the excitation of radial 
and non-radial acoustic pulsations within the context of asteroseismology. We have already imple¬ 
mented a deformable grid algorithm which we are now using to model oblate stars. The next step 
also include further development of CHORUS as a Large-Eddy Simulation (LES) code, with mod¬ 
eling of subgrid-scale (SGS) motions based either on explicit turbulence models or on the implicit 
numerical dissipation intrinsic to the SDM method. 
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Unstructured mesh consisting of all hexahedral elements. The full spherical shell 
mesh (a) is generated by fusing eight segments as illustrated in (b) using the GAM¬ 
BIT software [22]- At present, nodal points are uniformly distributed in the radial 

direction but GAMBIT also allows for variable radial resolution. 

Layout of a standard element for a 3rd order SDM in 2D . 

Strong scaling of the CHORUS code on Yellowstone for 18.9M (testl) and 70.8M 
(test2) DOFs. The latter test in particular achieves over 92% effeciency for 12k 

processors. 

Deviations from angular momentum conservation expressed in terms of angular ve¬ 


locity variations according to Eq. (33). Sfl x ,6Q y , and Sfl z are sampled every 5000 


iterations just before the correction procedure and normalized to the rotation rate 

of the coordinate system Do. 

(a) Kinetic energy density and ( b ) growth rate for the Jupiter benchmark. 

(a) Mach number and (6) e for the Jupiter benchmark, plotted at the sampling time 
T s . In (6), the initial condition is also plotted for comparison (dash-dotted line). 

The CHORUS and ASH curves are nearly indistinguishable. 

Components of the radial energy flux for the Jupiter benchmark integrated over 
horizontal surfaces at the sampling time T s for (a) CHORUS and (b) ASH. All 

values are normalized by the total flux F* imposed at the bottom boundary. 

Components of the radial energy flux as in Fig. [7] but at a much later time, t ~ 1800 

days, computed with the ASH code, showing an equilibrium solution. 

Mollweide projections of radial velocity V r for the Jupiter benchmark at the horizon¬ 
tal surface r = 0.95i? o , taken for (a) CHORUS and (b) ASH at the sampling time 
T s . Red and blue tones denote the upflow and downflow as indicated by the color bar. 
Shown are the power spectra of (a) radial velocity V r , (b ) meridional velocity Vg , 
and (c) zonal velocity V), as a function of spherical harmonic degree t for the Jupiter 

benchmark at r = 0.98 R 0 and t = T s . 

Mean flows in the Jupiter benchmark from CHORUS (top row) and ASH (bottom 
row). Shown are (a, e) meridional circulation expressed in terms of the stream 
function 4/, with red and blue denoting clockwise and counterclockwise circulations 
respectively, (6, /) specific entropy perturbation S , (c, g) temperature perturbation 
T , and (d, h) differential rotation, expressed in terms of the angular velocity Q. The 

top and bottom plots in each column share the same color bar. 

(a) Kinetic energy density and (b) growth rate for the solar benchmark. 

(a) Mach number and ( b) normalized entropy gradient e for the solar benchmark at 
the sampling time T s . As in Fig. |6(b)[ the initial profile is included in frame (6), 
where the CHORUS and ASH curves at T s are slightly steeper and indistinguishable 

from one another. 

Components of the normalized, horizontally-integrated radial energy flux as in Fig. 

[7] but for the solar benchmark. Each snapshot corresponds to the sampling time T s . 
Mollweide projections of the radial velocity V r in the solar benchmark at r = 0.95i? o 
and t = T s for (a) CHORUS and (b) ASH. Red and blue tones denote the upflow 

and downflow as indicated by the color bar. 

Power spectra of (a) radial velocity V r , (b) meridional velocity Vg, and (c) zonal 


velocity as in Fig. 10 but for the solar benchmark at r = 0.98 R a and t = T S . 
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17 Mean flows as in FigJTT] but for the solar benchmark. Top and bottom rows cor¬ 
respond to CHORUS and ASH respectively and all quantities are averaged over 
longitude and over a two-day time interval beginning at T s (spanning 2.2 rotation 
periods), (a, e) mass flux stream function 'k with red and blue denoting clockwise 
and counterclockwise circulations respectively, (b, /) specific entropy perturbation 
S , (c, g) temperature perturbation T , and (d, h) angular velocity profile. The top 
and bottom plots in each column share the same color bar. 
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